%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function builds proxies for all objects and groups
% them two by two
% 
% nobjsx   number of objects in the x-direction
% nobjsy   number of objects in the x-direction
% contours list of all the contours
% MatsS    list of discretization matrices corresponding to all
%          contours
% kh       Helmholtz parameter.
% acc      accuracy for builsing low rank approximations
% ppl      required number of points per wavelength on the outer
%          contours
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [contours_up,matsA_up,nobjsx_new,nobjsy_new,Dplt,Css] = factor_level(nobjsx,nobjsy,contours,matsA,kh,acc,ppl,ifig,aparams)
if nargin > 7
  plt = true;
  Dplt=cell(2,1);% contours to be plotted later because otherwise 'overplotted' by lated plotting
  Css=cell(2,1);
else
  plt = false;
end
if nargin > 8
  anim = true;
else
  anim = false;
end

nobjs       = nobjsx*nobjsy  ;
contours_up = cell(nobjs/2,1);
matsA_up    = cell(nobjs/2,1);

if plt
  figure(ifig);
end
if anim
  clf; hold all
  axis off
  axis(aparams.axis)
  axis equal
end

if plt
  for i=1:nobjs
    Ctmp = contours{i};
    plot(Ctmp(1,:) ,Ctmp(4,:) ,'k.')
  end
end
if anim
  fname = sprintf(aparams.fname,aparams.isave);
  saveas(gcf,fname)
  aparams.isave = aparams.isave + 1;
end

% choose the grouping direction
if nobjsx > nobjsy;
  sx   = 2;
  sy   = 1;
  step = 1;
  nobjsx_new = nobjsx/2;
  nobjsy_new = nobjsy;
else
  step = nobjsx;
  sx   = 1;
  sy   = 2;
  nobjsx_new = nobjsx;
  nobjsy_new = nobjsy/2;
end

inew = 1;
for i=1:sx:nobjsx
  for j=1:sy:nobjsy
    ind1 = (j-1)*nobjsx + i;
    % Proxy for Gamma1
    % fprintf(1,'Skeleton for Gamma 1        : %d \n')
    [D,h_D,C_in] = get_rect(contours,ind1,kh,ppl);
    if plt
	  Dplt{1,i}=D;
      rect = plot(D(1,:),D(2,:),'r','linewidth',2);
    end
    if ~isempty(C_in)
      if plt
	cin = plot(C_in(1,:),C_in(4,:),'go','linewidth',2);
      end
      D = [D,C_in([1,4],:)];
    end
    if anim
      fname = sprintf(aparams.fname,aparams.isave);
      saveas(gcf,fname)
      aparams.isave = aparams.isave + 1;
    end
    
    [Is1,Cs1,V1,matP1,D] = get_skeleton(kh,contours{ind1}, ...
					   matsA{ind1},acc,D,h_D);
    if plt
      plot(Cs1(1,:),Cs1(4,:),'r.','markersize',10)
	  Css{1,i}=Cs1;
    end
    
    if anim
      fname = sprintf(aparams.fname,aparams.isave);
      saveas(gcf,fname)
      aparams.isave = aparams.isave + 1;
      delete(rect)
    end
    
    if ~isempty(C_in) && anim
      delete(cin)
    end
    
    ind2 = ind1 + step;
    
    % fprintf(1,'Skeleton for Gamma 2        : %d \n')
    [D,h_D,C_in] = get_rect(contours,ind2,kh,ppl);
    if plt
	  Dplt{2,i}=D;
      rect = plot(D(1,:),D(2,:),'r','linewidth',2);
    end
    if ~isempty(C_in)
      if plt
	cin = plot(C_in(1,:),C_in(4,:),'go','linewidth',2);
      end
      D = [D,C_in([1,4],:)];
    end
    
    if anim
      fname = sprintf(aparams.fname,aparams.isave);
      saveas(gcf,fname)
      aparams.isave = aparams.isave + 1;
    end
    
    [Is2,Cs2,V2,matP2,D] = get_skeleton(kh,contours{ind2}, ...
					   matsA{ind2},acc,D,h_D);
    if plt
      plot(Cs2(1,:),Cs2(4,:),'r.','markersize',10)
      Css{2,i}=Cs2;
    end
    
    if anim
      fname = sprintf(aparams.fname,aparams.isave);
      saveas(gcf,fname)
      aparams.isave = aparams.isave + 1;
      delete(rect)
    end
    
    if ~isempty(C_in) && anim
      delete(cin)
    end
    
    % Use the skeletons
    [Cnew,matAnew] = merge_skeletons(kh,Cs1,matP1,Cs2,matP2);
    
    contours_up{inew} = Cnew;
    matsA_up{inew}    = matAnew;
    inew = inew +1;
    
    if anim
      plot(Cnew(1,:),Cnew(4,:),'.','markersize',2)
      fname = sprintf(aparams.fname,aparams.isave);
      saveas(gcf,fname)
      aparams.isave = aparams.isave + 1;
    end
  end
end

return